Seasonal shuffling
SEASONAL TEMPERATURES
Figure 1: Temperatures are seasonal
temp3=read.csv("/Users/Daisy/Documents/CRFL/Tissue vectoring/qPCR/R/Emerald_vector_cores/foweytemps5yrs.csv")
temp3$date=as.Date(with(temp3, paste('2019',MM, DD,sep="-")), "%Y-%m-%d")
temp3$WTMP=as.numeric(paste(temp3$WTMP))
temp=read.csv("/Users/Daisy/Documents/CRFL/Tissue vectoring/qPCR/R/Emerald_vector_cores/temp.csv")
temp$Date=as.Date(temp$Date,format = "%d/%m/%Y")
ggplot()+stat_summary(data=temp3,aes(x=date,y=WTMP,fill='darkgreen'),fun.data=mean_se,geom='ribbon',alpha=0.6)+theme_minimal(base_size = 15)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())+
labs(x='Month (2019-20)',y='Temperature/°C')+stat_summary(data=temp,aes(x=Date,y=Temp,fill='chartreuse2'),fun.data=mean_se,geom='ribbon',alpha=0.6)+scale_fill_identity(guide='legend',breaks=c('darkgreen','chartreuse2'),labels=c('Fowey Rock 2015-19 average','Emerald Reef 2019'))+guides(fill=guide_legend(title=''))+theme(legend.position = c(0.6,0.3))+scale_x_date(date_labels = '%b',limits=as.Date(c('2019-01-01','2019-12-31')))
#ggsave('temp3.pdf',device='pdf', width=7,height=5) #add labels to figure to show two collection timepointsINITIAL SYMBIONT COMMUNITIES
Figure 2a: M. cavernosa symbionts per host cell increased from April to October.
Add arrows to show that points above y=0 show higher in april, below line shows higher in April. Datapoints show the change in average S:H per colony. Stat_summary shows median and IQR for each species.
#calculate seasonal per colony change in symbiont density (ie the average symbiont density between the two cores taken from each colony per batch)
batches$pre_sh[which(batches$pre_sh==0)]=NA
Apr_Oct_sh= batches %>%
group_by(Colony,Batch) %>%
summarise_at(vars(pre_sh), funs(mean(., na.rm=TRUE)))
Apr_sh= filter(Apr_Oct_sh, Batch=='April')
Apr_sh$Apr_pre_sh=paste(Apr_sh$pre_sh)
Oct_sh= filter(Apr_Oct_sh, Batch=='October')
Oct_sh$Oct_pre_sh=paste(Oct_sh$pre_sh)
Apr_Oct_sh_change=join(Apr_sh,Oct_sh,by='Colony')
Apr_Oct_sh_change=Apr_Oct_sh_change[,c(1,4,7)]
Apr_Oct_sh_change$Apr_pre_sh=as.numeric(Apr_Oct_sh_change$Apr_pre_sh)
Apr_Oct_sh_change$Oct_pre_sh=as.numeric(Apr_Oct_sh_change$Oct_pre_sh)
Apr_Oct_sh_change$pre_sh_change= Apr_Oct_sh_change$Oct_pre_sh- Apr_Oct_sh_change$Apr_pre_sh
Apr_Oct_sh_change$rel_sh_change= ((Apr_Oct_sh_change$Oct_pre_sh- Apr_Oct_sh_change$Apr_pre_sh)/Apr_Oct_sh_change$Apr_pre_sh)*100
Apr_Oct_sh_change$octtoapr<- Apr_Oct_sh_change$Oct_pre_sh/Apr_Oct_sh_change$Apr_pre_sh
Apr_Oct_sh_change$Species=factor(Apr_Oct_sh_change$Colony)
Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
from=c('100','2','27','28','39','67','68','71','72','87'),to=(rep('M.cavernosa',times=10)))
Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
from=c('13','21','34','36','4','62','64','65','66','81'),to=(rep('O.faveolata',times=10)))
Apr_Oct_sh_change$Species= mapvalues(Apr_Oct_sh_change$Species,
from=c('16','18','19','22','23','26','3','35','41','48'),to=(rep('S.siderea',times=10)))
# this dataframe shows the average change in s:h per colony
##test statistical significance in seasonal change in S:H
hist(log10(batches$pre_sh)) #log10 transformed data then used linear mixed effects model on s:h data to test effect of batch within each species, with colony as a random factor
batches$transf_presh= log10(batches$pre_sh) # transform the response variable
mcavpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='M.cavernosa'))
rc_resids<- compute_redres(mcavpreshmod)
resids<- subset(batches,batches$Species=='M.cavernosa')
resids$logpresh<- log10(resids$pre_sh)
resids<-resids[,c(12)]
logpre<-resids[-c(5,6)]
resids<- data.frame(logpre, rc_resids)
plot_resqq(mcavpreshmod) # check residuals are normally distributed
mcavpreshmod2<- as_lmerModLmerTest(mcavpreshmod)
summary(mcavpreshmod2) #p for batch, p=0.001126** ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "M.cavernosa")
##
## REML criterion at convergence: 56
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.17638 -0.65717 -0.08901 0.55392 2.38020
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 3.240e-08 0.00018
## Residual 2.357e-01 0.48546
## Number of obs: 38, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.3177 0.1144 35.9720 -20.25 < 2e-16 ***
## BatchOctober 0.5583 0.1577 35.9823 3.54 0.00113 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.725# now let's get the batch parameter estimates and CIs:
mcavemm.sh<- emmeans(mcavpreshmod, specs=revpairwise~Batch)
summary(mcavemm.sh)## $emmeans
## Batch emmean SE df lower.CL upper.CL
## April -2.32 0.114 36 -2.55 -2.09
## October -1.76 0.109 36 -1.98 -1.54
##
## Degrees-of-freedom method: satterthwaite
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## October - April 0.558 0.158 36 3.540 0.0011
##
## Note: contrasts are still on the log10 scale
## Degrees-of-freedom method: satterthwaitemcavemmsh_contrasts<- mcavemm.sh$contrasts %>%
confint()%>%
as.data.frame() # NB these results are given on a log10 scale
#also include April prediction in this dataframe for later calculations
mcavemmsh_contrasts<-as.data.frame(c(mcavemmsh_contrasts,mcavemm.sh$emmeans[1]))
ofavpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='O.faveolata'))
rc_resids<- compute_redres(ofavpreshmod)
resids<- subset(batches,batches$Species=='O.faveolata')
resids$logpresh<- log10(resids$pre_sh)
resids<-resids[,c(12)]
logpre<-resids
resids<- data.frame(logpre, rc_resids)
plot_resqq(mcavpreshmod) # check residuals are normally distributed
ofavpreshmod2<- as_lmerModLmerTest(ofavpreshmod)
summary(ofavpreshmod2) #p for batch, p=0.938## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "O.faveolata")
##
## REML criterion at convergence: 29.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3754 -0.5262 0.1305 0.6129 1.6623
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.12853 0.3585
## Residual 0.07002 0.2646
## Number of obs: 35, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.491446 0.128795 11.205854 -11.580 1.4e-07 ***
## BatchOctober 0.007186 0.091378 24.399854 0.079 0.938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.320# now let's get the batch parameter estimates and CIs:
ofavemm.sh<- emmeans(ofavpreshmod, specs=revpairwise~Batch)
summary(ofavemm.sh)## $emmeans
## Batch emmean SE df lower.CL upper.CL
## April -1.49 0.129 11.2 -1.77 -1.21
## October -1.48 0.132 12.2 -1.77 -1.20
##
## Degrees-of-freedom method: satterthwaite
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## October - April 0.00719 0.0914 24.4 0.079 0.9380
##
## Note: contrasts are still on the log10 scale
## Degrees-of-freedom method: satterthwaiteofavemmsh_contrasts<- ofavemm.sh$contrasts %>%
confint()%>%
as.data.frame() # NB these results are given on a log10 scale
#also include April prediction in this dataframe for later calculations
ofavemmsh_contrasts<-as.data.frame(c(ofavemmsh_contrasts,ofavemm.sh$emmeans[1]))
ssidpreshmod=lmer(log10(pre_sh)~Batch+(1|Colony),data=filter(batches,batches$Species=='S.siderea'))
plot(ssidpreshmod)
rc_resids<- compute_redres(ssidpreshmod)
resids<- subset(batches,batches$Species=='S.siderea')
resids$logpresh<- log10(resids$pre_sh)
resids<-resids[,c(12)]
logpre<-resids[-c(19,20,39,40)]
resids<- data.frame(logpre, rc_resids)
plot_resqq(mcavpreshmod) # check residuals are normally distributed
ssidpreshmod2<- as_lmerModLmerTest(ssidpreshmod)
summary(ssidpreshmod2) #p for batch, p=0.04833*## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(pre_sh) ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "S.siderea")
##
## REML criterion at convergence: 55
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.53246 -0.34681 0.00358 0.49889 1.82110
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.2705 0.5201
## Residual 0.1524 0.3903
## Number of obs: 36, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.0287 0.1963 10.0494 -15.431 2.51e-08 ***
## BatchOctober -0.2696 0.1301 26.0000 -2.072 0.0483 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.331# now let's get the batch parameter estimates and CIs:
ssidemm.sh<- emmeans(ssidpreshmod, specs=revpairwise~Batch)
summary(ssidemm.sh)## $emmeans
## Batch emmean SE df lower.CL upper.CL
## April -3.03 0.196 10.1 -3.47 -2.59
## October -3.30 0.196 10.1 -3.74 -2.86
##
## Degrees-of-freedom method: satterthwaite
## Results are given on the log10 (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## October - April -0.27 0.13 26 -2.072 0.0483
##
## Note: contrasts are still on the log10 scale
## Degrees-of-freedom method: satterthwaitessidemmsh_contrasts<- ssidemm.sh$contrasts %>%
confint()%>%
as.data.frame() # NB these results are given on a log10 scale
#also include April prediction in this dataframe for later calculations
ssidemmsh_contrasts<-as.data.frame(c(ssidemmsh_contrasts,ssidemm.sh$emmeans[1]))
#collate the dataframe
sh_emmcontrasts<-rbind(mcavemmsh_contrasts,ofavemmsh_contrasts,ssidemmsh_contrasts)
sh_emmcontrasts$Species<-c('M.cavernosa','O.faveolata','S.siderea')
sh_emmcontrasts$Species<-as.factor(sh_emmcontrasts$Species)
sh_emmcontrasts$change_untransformed= 10^(sh_emmcontrasts$estimate)
sh_emmcontrasts$lowerCI_change_untransformed= 10^(sh_emmcontrasts$lower.CL)
sh_emmcontrasts$upperCI_change_untransformed= 10^(sh_emmcontrasts$upper.CL)
sh_emmcontrasts$April_untransformed= 10^(sh_emmcontrasts$emmean)
#plot set-up
ofav=expression(paste(italic("O. faveolata")))
ssid=expression(paste(italic("S. siderea")))
mcav=expression(paste(italic("M. cavernosa")))
#PLOT 2A.SEASONAL CHANGE IN INITIAL S:H
ggplot()+
geom_point(data=sh_emmcontrasts, aes(x=Species,y=(change_untransformed), colour=Species), size=2)+
geom_point(data=Apr_Oct_sh_change, aes(x=Species,y=octtoapr, colour=Species), size=0.5, position=position_jitter(width=0.1))+
geom_errorbar(data=sh_emmcontrasts, aes(x=Species, ymin=(lowerCI_change_untransformed), ymax=(upperCI_change_untransformed), colour=Species), size=0.5, width=0.2)+
theme_minimal(base_size = 15)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())+
labs(y='October : April symbionts per host cell', x='')+
scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(add=c(0.4,1.2)))+
scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
guides(colour=F)+
geom_hline(yintercept=1,linetype='dashed')+
scale_y_continuous(breaks=c(0.5,1,2,4,6,8))+
coord_cartesian(ylim=c(0.3,8))
#ggsave('initialSH_seasonal_emmeans.pdf',device='pdf',width=7,height=5) # add 'higher in Oct'/'higher in Apr' labels to figure
#UPDATES FOR ROSS: This is now the emmeans predicted seasonal differnece, with 95% confidence intervals. Note some individual mcav datapoints are beyond the limits of the graph but accounted for in the predicted means and statistical calculations. Figure 2b: Proportion Durusdinium did not change between April and October for any of the 3 coral species.
batch_comparevi=
ggplot()+
theme_minimal(base_size = 15)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())+
labs(y='Initial proportion *Durusdinium*')+
scale_fill_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
geom_violin(data=batches,aes(x=Species,y=pre_propD,colour=Batch),scale = 'width')+
geom_dotplot(data=batches,bins=30,binaxis='y',dotsize=0.7,stackratio=0.5,stackdir='center',stackgroups=F, position='dodge',aes(x= Species,y=pre_propD,fill=Batch))+
scale_colour_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
scale_x_discrete(labels=c(mcav,ofav,ssid))+
theme(axis.title.y.left = element_markdown(), legend.position = c(0.2,0.5))+
labs(x='')+
theme(legend.title = element_text(size=0))
batch_comparevi
#ggsave('batchcompare_propd_vi_dot.pdf',device='pdf',width=7,height=5)
hist(batches$pre_propD)
mcavpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='M.cavernosa'), family = 'binomial')
summary(mcavpredmod) # p=1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "M.cavernosa")
##
## AIC BIC logLik deviance df.resid
## 6.0 11.1 0.0 0.0 37
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## 0 0 0 0 2904892
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0 0
## Number of obs: 40, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.806e+01 1.501e+07 0 1
## BatchOctober -3.033e+02 2.122e+07 0 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.707
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularofavpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='O.faveolata'), family = 'binomial')
summary(ofavpredmod) # p=0.2## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "O.faveolata")
##
## AIC BIC logLik deviance df.resid
## 31.9 36.5 -12.9 25.9 32
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.76647 -0.04036 -0.00779 0.25178 1.00309
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 80.35 8.964
## Number of obs: 35, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.153 3.793 -1.622 0.105
## BatchOctober -3.286 2.769 -1.187 0.235
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr 0.343ssidpredmod=glmer(pre_propD~Batch+(1|Colony),data=filter(batches,batches$Species=='S.siderea'), family = 'binomial')
summary(ssidpredmod)# p0.09## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pre_propD ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "S.siderea")
##
## AIC BIC logLik deviance df.resid
## 36.2 41.0 -15.1 30.2 33
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9098 -0.3598 -0.0362 0.2155 0.9639
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 24.15 4.914
## Number of obs: 36, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.039 2.132 -0.487 0.6259
## BatchOctober 4.593 2.713 1.693 0.0905 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.425#UPDATES FOR ROSS: Make it clear that this plot is showing raw data, not predictive model. glmer is now 'family quasibinomial'. BLEACHING SENSITIVITY & RATE
Figure 3: M. cavernosa showed higher sensitivity to bleaching in October compared to April.
blch_sensitivity=blch_sensitivity[-c(31,32,48,59,56,57,75),] #remove 're'_drop_sh' NAs, and 5 cores that increased symbiont density
mcavsensitivity=subset(blch_sensitivity,Species=='M.cavernosa')
mcavsensitivity$transformedshdrop<- (mcavsensitivity$rel_drop_sh)^2
plot(mcavsensitivity$rel_drop_sh~mcavsensitivity$rel_drop_y2)
mcavblchresmod=glmer((rel_drop_sh^2)~rel_drop_y2*Batch+(1|Colony),data=mcavsensitivity)
plot_resqq(mcavblchresmod)
summary(mcavblchresmod)## Linear mixed model fit by REML ['lmerMod']
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + (1 | Colony)
## Data: mcavsensitivity
##
## REML criterion at convergence: 363.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.90415 -0.22251 0.04061 0.26312 2.78896
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 30138 173.6
## Residual 746489 864.0
## Number of obs: 25, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -5141.38 2432.88 -2.113
## rel_drop_y2 -97.15 34.76 -2.795
## Batch2 15479.84 2674.24 5.789
## rel_drop_y2:Batch2 106.79 40.53 2.635
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2
## rel_drop_y2 0.993
## Batch2 -0.912 -0.906
## rl_drp_2:B2 -0.856 -0.862 0.985mcavblchresmod2<- as_lmerModLmerTest(mcavblchresmod)
summary(mcavblchresmod2)#the interaction between batch and drop in y2 is significant (p=0.0156).## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + (1 | Colony)
## Data: mcavsensitivity
##
## REML criterion at convergence: 363.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.90415 -0.22251 0.04061 0.26312 2.78896
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 30138 173.6
## Residual 746489 864.0
## Number of obs: 25, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -5141.38 2432.88 21.00 -2.113 0.0467 *
## rel_drop_y2 -97.15 34.76 21.00 -2.795 0.0109 *
## Batch2 15479.84 2674.24 20.94 5.789 9.67e-06 ***
## rel_drop_y2:Batch2 106.79 40.53 20.73 2.635 0.0156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2
## rel_drop_y2 0.993
## Batch2 -0.912 -0.906
## rl_drp_2:B2 -0.856 -0.862 0.985ofavsensitivity<-subset(blch_sensitivity,Species=='O.faveolata')
plot(ofavsensitivity$rel_drop_sh~ofavsensitivity$rel_drop_y2)
ofavblchresmod=glmer((rel_drop_sh^2)~rel_drop_y2*Batch+InitialDom+(1|Colony),data=ofavsensitivity)
plot_resqq(ofavblchresmod)
summary(ofavblchresmod)## Linear mixed model fit by REML ['lmerMod']
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + InitialDom + (1 | Colony)
## Data: ofavsensitivity
##
## REML criterion at convergence: 335.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.48257 -0.28194 0.08217 0.54525 1.33783
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 195290 441.9
## Residual 756932 870.0
## Number of obs: 24, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 7700.76 1267.88 6.074
## rel_drop_y2 -21.06 22.82 -0.923
## Batch2 3731.23 2596.95 1.437
## InitialDomnond 469.74 584.43 0.804
## rel_drop_y2:Batch2 64.05 41.25 1.553
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2 IntlDm
## rel_drop_y2 0.926
## Batch2 -0.576 -0.704
## InitilDmnnd 0.036 0.349 -0.525
## rl_drp_2:B2 -0.593 -0.742 0.988 -0.539ofavblchresmod2<- as_lmerModLmerTest(ofavblchresmod)
summary(ofavblchresmod2) #when separated by symbiont genus, the interaction between bacth and drop in y2 is insignificant (p=0.137)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + InitialDom + (1 | Colony)
## Data: ofavsensitivity
##
## REML criterion at convergence: 335.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.48257 -0.28194 0.08217 0.54525 1.33783
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 195290 441.9
## Residual 756932 870.0
## Number of obs: 24, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7700.76 1267.88 18.99 6.074 7.7e-06 ***
## rel_drop_y2 -21.06 22.82 18.40 -0.923 0.368
## Batch2 3731.23 2596.95 18.97 1.437 0.167
## InitialDomnond 469.74 584.43 11.22 0.804 0.438
## rel_drop_y2:Batch2 64.05 41.25 19.00 1.553 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2 IntlDm
## rel_drop_y2 0.926
## Batch2 -0.576 -0.704
## InitilDmnnd 0.036 0.349 -0.525
## rl_drp_2:B2 -0.593 -0.742 0.988 -0.539ssidsensitivity<-subset(blch_sensitivity,Species=='S.siderea')
plot(ssidsensitivity$rel_drop_sh~ssidsensitivity$rel_drop_y2)
ssidblchresmod=glmer((rel_drop_sh^2)~rel_drop_y2*Batch+InitialDom+(1|Colony),data=ssidsensitivity)
plot_resqq(ssidblchresmod)
summary(ssidblchresmod)## Linear mixed model fit by REML ['lmerMod']
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + InitialDom + (1 | Colony)
## Data: ssidsensitivity
##
## REML criterion at convergence: 269.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7786 -0.3998 0.2009 0.6135 1.1655
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0 0.0
## Residual 868824 932.1
## Number of obs: 20, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8273.79 1130.84 7.316
## rel_drop_y2 -15.86 22.12 -0.717
## Batch2 1418.57 1413.00 1.004
## InitialDomnond 218.95 459.78 0.476
## rel_drop_y2:Batch2 21.21 27.17 0.781
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2 IntlDm
## rel_drop_y2 0.942
## Batch2 -0.747 -0.757
## InitilDmnnd -0.224 0.014 -0.060
## rl_drp_2:B2 -0.733 -0.816 0.947 -0.163
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularssidblchresmod2<- as_lmerModLmerTest(ssidblchresmod)
summary(ssidblchresmod2) #when separated by symbiont genus, the interaction between bacth and drop in y2 is insignificant (p=0.447)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (rel_drop_sh^2) ~ rel_drop_y2 * Batch + InitialDom + (1 | Colony)
## Data: ssidsensitivity
##
## REML criterion at convergence: 269.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7786 -0.3998 0.2009 0.6135 1.1655
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0 0.0
## Residual 868824 932.1
## Number of obs: 20, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8273.79 1130.84 15.00 7.316 2.54e-06 ***
## rel_drop_y2 -15.86 22.12 15.00 -0.717 0.484
## Batch2 1418.57 1413.00 15.00 1.004 0.331
## InitialDomnond 218.95 459.78 15.00 0.476 0.641
## rel_drop_y2:Batch2 21.21 27.17 15.00 0.781 0.447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rl_d_2 Batch2 IntlDm
## rel_drop_y2 0.942
## Batch2 -0.747 -0.757
## InitilDmnnd -0.224 0.014 -0.060
## rl_drp_2:B2 -0.733 -0.816 0.947 -0.163
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularmcavsensitivity$predicted<- predict(mcavblchresmod)
mcavsensitivity$residuals<- residuals(mcavblchresmod)
mcavsensitivity$untransformed_predicted<- -(mcavsensitivity$predicted^0.5) #all changes are negative, reverse transform sqaured response variable.
ggplot()+
geom_point(data=mcavsensitivity,aes(x=rel_drop_y2,y=rel_drop_sh,colour=Batch))+
geom_smooth(data=mcavsensitivity,method='lm',aes(x=rel_drop_y2,y=untransformed_predicted,color=Batch),show.legend=F,se=T, alpha=0.5, size=0.5)+
coord_cartesian(clip='off', ylim=c(-100,0))+
theme_minimal(base_size = 13)%+replace%
theme(panel.border = element_rect(colour='grey20',size=0.5,fill=NA))+
scale_colour_manual(values=c('#3C5488B2','#00A087B2'),labels=c('April','October'))+
theme(panel.spacing = unit(1.2, "lines"),legend.position = 'right')+
theme(legend.title=element_text(size=0))+
theme(legend.position = c(0.8,0.5))+
labs(x='% change in Fv/Fm',y='% change in symbionts per host cell')
#ggsave('bleaching_sensitivity_batches.pdf',device='pdf',width=7,height=5)
#UPDATES FOR ROSS: This model is now also a mixed effects model, performed on transformed response data. A glmm was used instead of a lmm due to uneven sample sizes of from each colony in each treatemnt due to exclusion of control cores. This now shows the linear regression and 95% confidence interval for predicted values from the model. The plot shows MCAV only but stats have been done for all three species. Those with >0.75 initial prop D were grouped and those with <0.25 prop D were grouped (this excluded very few datapoints). Since cores were pulled and sampled after exposure to differnet DHW, there is variation in DHW exposure, as displayed on the x-axis below.
Figure 4: Coral species and dominant symbiont genus both shape the trajectories of symbiont loss and recovery over bleaching.
allbleach=rbind(ofav_DHW,ssid_DHW,mcav_DHW)
allbleach$Species=factor(allbleach$Species,levels=c('M.cavernosa','O.faveolata','S.siderea'))
allrecov=rbind(mcav_recov, ofav_recov, ssid_recov)
allrecov$Species=factor(allrecov$Species,levels=c('M.cavernosa','O.faveolata','S.siderea'))
allrecov$InitialDom=revalue(allrecov$InitialDom,c('D'='d','B'='nond','C'='nond'))
allrecov$InitialDom=factor(allrecov$InitialDom)
allbleach$InitialDom=revalue(allbleach$InitialDom,c('D'='d','B'='nond','C'='nond'))
allbleach$InitialDom=factor(allbleach$InitialDom)
allbleach$Batch=factor(allbleach$Batch)
allrecov$Batch=factor(allrecov$Batch)
allbleach$Colony=factor(allbleach$Colony)
allrecov$Colony=factor(allrecov$Colony)
allbleach<-allbleach[-c(167),]
allrecov<-allrecov[-c(86,239),]#model highlighted this one datapoint as an outlier in both dataframes
bleachingmod<-glm(Shprop~DHW*Species*InitialDom, data=allbleach, family='quasipoisson') #plot data with a quassipoisson model
bleachingfit= with(summary(bleachingmod), 1 - deviance/null.deviance)
bleachingfit #0.8533836, R2 for plotted quasipoisson model with batch, 0.6939 without batch## [1] 0.6938675bleachingmod2<-glmer(Shprop~DHW*Species*InitialDom+(1|Colony), data=allbleach)
plot_resqq(bleachingmod2) #perform statistical tests on mixed effects model
bleachingmod3<-as_lmerModLmerTest(bleachingmod2)
summary(bleachingmod3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * Species * InitialDom + (1 | Colony)
## Data: allbleach
##
## REML criterion at convergence: 62.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1750 -0.4013 -0.0358 0.3297 5.1995
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.006468 0.08042
## Residual 0.060175 0.24531
## Number of obs: 147, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.054054 0.125491 77.511446 8.399
## DHW -0.059963 0.017904 120.365702 -3.349
## SpeciesO.faveolata -0.057977 0.159516 67.688528 -0.363
## SpeciesS.siderea -0.099000 0.102079 67.803938 -0.970
## InitialDomnond -0.023200 0.113848 86.620959 -0.204
## DHW:SpeciesO.faveolata -0.009859 0.020228 119.185291 -0.487
## DHW:SpeciesS.siderea 0.002114 0.016147 121.484396 0.131
## DHW:InitialDomnond -0.044485 0.015250 121.587128 -2.917
## SpeciesO.faveolata:InitialDomnond -0.033899 0.164764 68.629139 -0.206
## DHW:SpeciesO.faveolata:InitialDomnond 0.023099 0.019905 120.391812 1.161
## Pr(>|t|)
## (Intercept) 1.65e-12 ***
## DHW 0.00108 **
## SpeciesO.faveolata 0.71740
## SpeciesS.siderea 0.33557
## InitialDomnond 0.83901
## DHW:SpeciesO.faveolata 0.62690
## DHW:SpeciesS.siderea 0.89605
## DHW:InitialDomnond 0.00421 **
## SpeciesO.faveolata:InitialDomnond 0.83760
## DHW:SpeciesO.faveolata:InitialDomnond 0.24814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DHW SpcsO. SpcsS. IntlDm DHW:SpO. DHW:SS DHW:ID SO.:ID
## DHW -0.592
## SpecsO.fvlt -0.787 0.466
## SpecisS.sdr -0.798 0.536 0.627
## InitilDmnnd -0.907 0.507 0.714 0.639
## DHW:SpcsO.f 0.524 -0.885 -0.590 -0.475 -0.449
## DHW:SpcsS.s 0.482 -0.901 -0.379 -0.597 -0.370 0.798
## DHW:IntlDmn 0.540 -0.852 -0.425 -0.440 -0.596 0.754 0.701
## SpcsO.fv:ID 0.627 -0.351 -0.862 -0.442 -0.691 0.483 0.256 0.412
## DHW:SpO.:ID -0.414 0.653 0.507 0.337 0.456 -0.798 -0.537 -0.766 -0.597
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients anova(bleachingmod3)#no significant difference in bleaching (DHW interaction with species) between species p=0.818834, but a significant interaction between DHW and initial symbiont type (p=0.00421), with nond-hosting corals bleaching more. ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## DHW 22.5923 22.5923 1 118.404 375.4414 < 2.2e-16 ***
## Species 0.0952 0.0476 2 56.230 0.7907 0.458514
## InitialDom 0.0143 0.0143 1 68.629 0.2375 0.627558
## DHW:Species 0.0241 0.0120 2 120.106 0.2002 0.818834
## DHW:InitialDom 0.6590 0.6590 1 120.392 10.9518 0.001234 **
## Species:InitialDom 0.0025 0.0025 1 68.629 0.0423 0.837599
## DHW:Species:InitialDom 0.0810 0.0810 1 120.392 1.3468 0.248139
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#now looking at any batch differences within each coral species:
mcavbleachingmod<-glmer(Shprop~DHW*Batch+(1|Colony), data=filter(allbleach,allbleach$Species=='M.cavernosa'))
plot_resqq(mcavbleachingmod)
mcavbleachingmod2<-as_lmerModLmerTest(mcavbleachingmod)
summary(mcavbleachingmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * Batch + (1 | Colony)
## Data: filter(allbleach, allbleach$Species == "M.cavernosa")
##
## REML criterion at convergence: -5.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4227 -0.0828 0.0082 0.0650 5.2633
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.00000 0.0000
## Residual 0.03716 0.1928
## Number of obs: 55, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.01595 0.05271 51.00000 19.274 < 2e-16 ***
## DHW -0.04102 0.01307 51.00000 -3.138 0.00283 **
## Batch2 -0.02849 0.07227 51.00000 -0.394 0.69507
## DHW:Batch2 -0.08481 0.01590 51.00000 -5.335 2.21e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DHW Batch2
## DHW -0.682
## Batch2 -0.729 0.497
## DHW:Batch2 0.561 -0.822 -0.682
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular anova(mcavbleachingmod2) #significant interaction between DHW and batch on proportion of symbionts retained, p=2.21e-06, with the october batch losing more.## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## DHW 4.0937 4.0937 1 51 110.1738 2.410e-14 ***
## Batch 0.0058 0.0058 1 51 0.1554 0.6951
## DHW:Batch 1.0577 1.0577 1 51 28.4652 2.207e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ofavbleachingmod<-glmer(Shprop~DHW*InitialDom*Batch+(1|Colony), data=filter(allbleach,allbleach$Species=='O.faveolata'))
plot_resqq(ofavbleachingmod)
ofavbleachingmod2<-as_lmerModLmerTest(ofavbleachingmod)
summary(ofavbleachingmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * InitialDom * Batch + (1 | Colony)
## Data: filter(allbleach, allbleach$Species == "O.faveolata")
##
## REML criterion at convergence: -2.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.14099 -0.30534 0.05302 0.62284 1.61004
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.008981 0.09477
## Residual 0.019452 0.13947
## Number of obs: 48, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.989014 0.088909 18.437381 11.124 1.29e-09 ***
## DHW -0.072441 0.007740 30.462625 -9.360 1.80e-10 ***
## InitialDomnond -0.038492 0.105732 17.426880 -0.364 0.7202
## Batch2 0.021250 0.099653 31.174272 0.213 0.8325
## DHW:InitialDomnond -0.025863 0.010215 31.202369 -2.532 0.0166 *
## DHW:Batch2 0.004669 0.010722 30.423269 0.435 0.6663
## InitialDomnond:Batch2 -0.025393 0.125355 33.537940 -0.203 0.8407
## DHW:InitialDomnond:Batch2 0.006091 0.014644 30.797769 0.416 0.6803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DHW IntlDm Batch2 DHW:InD DHW:B2 InD:B2
## DHW -0.548
## InitilDmnnd -0.841 0.460
## Batch2 -0.557 0.483 0.468
## DHW:IntlDmn 0.415 -0.758 -0.529 -0.366
## DHW:Batch2 0.395 -0.722 -0.332 -0.691 0.547
## IntlDmnn:B2 0.443 -0.384 -0.535 -0.795 0.454 0.550
## DHW:IntD:B2 -0.289 0.529 0.366 0.506 -0.695 -0.732 -0.675 anova(ofavbleachingmod2) #no significant batchxDHW interaction, p=0.6663 ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## DHW 9.5879 9.5879 1 31.033 492.9089 < 2.2e-16 ***
## InitialDom 0.0064 0.0064 1 9.848 0.3268 0.580365
## Batch 0.0004 0.0004 1 33.538 0.0186 0.892267
## DHW:InitialDom 0.1878 0.1878 1 31.033 9.6552 0.004019 **
## DHW:Batch 0.0216 0.0216 1 30.798 1.1102 0.300241
## InitialDom:Batch 0.0008 0.0008 1 33.538 0.0410 0.840695
## DHW:InitialDom:Batch 0.0034 0.0034 1 30.798 0.1730 0.680328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ssidbleachingmod<-glmer(Shprop~DHW*InitialDom*Batch+(1|Colony), data=filter(allbleach,allbleach$Species=='S.siderea'))
plot_resqq(ssidbleachingmod)
ssidbleachingmod2<-as_lmerModLmerTest(ssidbleachingmod)
summary(ssidbleachingmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ DHW * InitialDom * Batch + (1 | Colony)
## Data: filter(allbleach, allbleach$Species == "S.siderea")
##
## REML criterion at convergence: 46.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9957 -0.6123 0.1110 0.4286 2.9641
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02324 0.1524
## Residual 0.07560 0.2750
## Number of obs: 44, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.9753174 0.1423833 29.2629965 6.850 1.52e-07 ***
## DHW -0.0713355 0.0133116 30.1054814 -5.359 8.37e-06 ***
## InitialDomnond -0.0182845 0.1833807 31.6026207 -0.100 0.9212
## Batch2 -0.0417829 0.1596818 32.1581895 -0.262 0.7952
## DHW:InitialDomnond -0.0397296 0.0219691 30.8016710 -1.808 0.0803 .
## DHW:Batch2 0.0229258 0.0176451 30.3789459 1.299 0.2036
## InitialDomnond:Batch2 -0.0238694 0.2719345 32.5330885 -0.088 0.9306
## DHW:InitialDomnond:Batch2 0.0002652 0.0370123 30.3510824 0.007 0.9943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) DHW IntlDm Batch2 DHW:InD DHW:B2 InD:B2
## DHW -0.601
## InitilDmnnd -0.756 0.466
## Batch2 -0.713 0.534 0.578
## DHW:IntlDmn 0.359 -0.606 -0.566 -0.330
## DHW:Batch2 0.464 -0.756 -0.365 -0.679 0.459
## IntlDmnn:B2 0.411 -0.313 -0.521 -0.596 0.407 0.401
## DHW:IntD:B2 -0.218 0.360 0.339 0.328 -0.586 -0.477 -0.667 anova(ssidbleachingmod2)#no significant effect of batch on the relatinship between DHW an symbiont loss (interaction), p=0.2036## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## DHW 5.5006 5.5006 1 30.940 72.7582 1.256e-09 ***
## InitialDom 0.0026 0.0026 1 25.563 0.0349 0.85320
## Batch 0.0120 0.0120 1 32.269 0.1593 0.69239
## DHW:InitialDom 0.3402 0.3402 1 30.915 4.4998 0.04202 *
## DHW:Batch 0.1175 0.1175 1 30.334 1.5546 0.22200
## InitialDom:Batch 0.0006 0.0006 1 32.533 0.0077 0.93059
## DHW:InitialDom:Batch 0.0000 0.0000 1 30.351 0.0001 0.99433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1recovmod<-glm(Shprop~Recov.Days*Species*InitialDom, data=allrecov, family='quasipoisson')
plot(recovmod)
recoveryfit= with(summary(recovmod), 1 - deviance/null.deviance)
recoveryfit #0.5469, R^2 for plotted quasipoisson model (without batch)## [1] 0.5469235recovmod2<-glmer(Shprop~Recov.Days*Species*InitialDom+(1|Colony), data=allrecov)
plot_resqq(recovmod2)
recovmod3<-as_lmerModLmerTest(recovmod2)
summary(recovmod3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * Species * InitialDom + (1 | Colony)
## Data: allrecov
##
## REML criterion at convergence: 1250.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7901 -0.2597 -0.0424 0.0247 10.8057
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.1652 0.4064
## Residual 17.4242 4.1742
## Number of obs: 217, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 3.67046 1.79581 120.33950
## Recov.Days -0.12507 0.03986 188.75632
## SpeciesO.faveolata -3.00438 2.22969 110.98030
## SpeciesS.siderea -3.44919 1.48290 124.40494
## InitialDomnond -3.65091 1.65218 130.04791
## Recov.Days:SpeciesO.faveolata 0.14223 0.05070 188.65474
## Recov.Days:SpeciesS.siderea 0.15282 0.03292 187.61232
## Recov.Days:InitialDomnond 0.16279 0.03721 188.75877
## SpeciesO.faveolata:InitialDomnond 2.96527 2.31684 114.40581
## Recov.Days:SpeciesO.faveolata:InitialDomnond -0.17354 0.05372 188.60438
## t value Pr(>|t|)
## (Intercept) 2.044 0.04315 *
## Recov.Days -3.138 0.00197 **
## SpeciesO.faveolata -1.347 0.18058
## SpeciesS.siderea -2.326 0.02164 *
## InitialDomnond -2.210 0.02887 *
## Recov.Days:SpeciesO.faveolata 2.805 0.00556 **
## Recov.Days:SpeciesS.siderea 4.642 6.48e-06 ***
## Recov.Days:InitialDomnond 4.376 2.00e-05 ***
## SpeciesO.faveolata:InitialDomnond 1.280 0.20318
## Recov.Days:SpeciesO.faveolata:InitialDomnond -3.231 0.00146 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Rcv.Dy SpcsO. SpcsS. IntlDm Rc.D:SO. R.D:SS R.D:ID SO.:ID
## Recov.Days -0.748
## SpecsO.fvlt -0.805 0.603
## SpecisS.sdr -0.825 0.625 0.664
## InitilDmnnd -0.920 0.702 0.741 0.694
## Rcv.Dys:SO. 0.588 -0.786 -0.746 -0.491 -0.552
## Rcv.Dys:SS. 0.624 -0.826 -0.503 -0.756 -0.543 0.649
## Rcv.Dys:InD 0.692 -0.933 -0.557 -0.536 -0.752 0.734 0.718
## SpcsO.fv:ID 0.656 -0.500 -0.867 -0.495 -0.713 0.655 0.387 0.536
## Rc.D:SO.:ID -0.479 0.646 0.643 0.371 0.521 -0.869 -0.497 -0.693 -0.746
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficientsanova(recovmod3)# symbiont recovery was significantly differnet between species (recovery days * symbiont density interaction), with ofav p=0.00556 and ssid p=6.48e-06 recovering more than mcav for a given amount of recovery (likely indicative of the delay in symbiont recovery seen in mcav). There was also a significant interaction between initial symbiont genus and recovery days p=2.00e-05, with corals initially not hosting Durusdinium recovering more with a given amount of recovery. ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Recov.Days 336.61 336.61 1 188.77 19.3184 1.844e-05 ***
## Species 81.15 40.57 2 103.20 2.3286 0.1025280
## InitialDom 61.05 61.05 1 114.41 3.5035 0.0637941 .
## Recov.Days:Species 304.78 152.39 2 188.33 8.7458 0.0002333 ***
## Recov.Days:InitialDom 139.59 139.59 1 188.60 8.0111 0.0051531 **
## Species:InitialDom 28.54 28.54 1 114.41 1.6381 0.2031778
## Recov.Days:Species:InitialDom 181.85 181.85 1 188.60 10.4364 0.0014578 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 mcavrecovmod<-glmer(Shprop~Recov.Days*Batch+(1|Colony), data=filter(allrecov,allrecov$Species=='M.cavernosa'))
plot_resqq(mcavrecovmod)
mcavrecovmod2<-as_lmerModLmerTest(mcavrecovmod)
summary(mcavrecovmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * Batch + (1 | Colony)
## Data: filter(allrecov, allrecov$Species == "M.cavernosa")
##
## REML criterion at convergence: 416.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5973 -0.4961 -0.1152 0.1772 5.0981
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.1195 0.3457
## Residual 8.3119 2.8830
## Number of obs: 82, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.44603 0.75590 66.15752 -0.590 0.55716
## Recov.Days 0.08342 0.01772 70.99172 4.706 1.21e-05 ***
## Batch2 0.23704 0.98147 74.25824 0.242 0.80982
## Recov.Days:Batch2 -0.06258 0.02139 70.90131 -2.926 0.00461 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Rcv.Dy Batch2
## Recov.Days -0.769
## Batch2 -0.753 0.592
## Rcv.Dys:Bt2 0.637 -0.829 -0.755 anova(mcavrecovmod2) #there was a significant interaction p=0.004611 of batch on the relationchip between symbiont density and recovery days, with the October batch recovering less. ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Recov.Days 197.495 197.495 1 70.951 23.7606 6.437e-06 ***
## Batch 0.485 0.485 1 74.258 0.0583 0.809820
## Recov.Days:Batch 71.158 71.158 1 70.901 8.5610 0.004611 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ofavrecovmod<-glmer(Shprop~Recov.Days*InitialDom*Batch+(1|Colony), data=filter(allrecov,allrecov$Species=='O.faveolata'))
plot_resqq(ofavrecovmod)
ofavrecovmod2<-as_lmerModLmerTest(ofavrecovmod)
summary(ofavrecovmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * InitialDom * Batch + (1 | Colony)
## Data: filter(allrecov, allrecov$Species == "O.faveolata")
##
## REML criterion at convergence: 245.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5557 -0.3442 -0.1182 0.0422 6.6435
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.000 0.000
## Residual 1.413 1.189
## Number of obs: 70, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.59242 0.56690 62.00000 1.045 0.300
## Recov.Days 0.02810 0.01714 62.00000 1.640 0.106
## InitialDomnond -0.63370 0.66320 62.00000 -0.956 0.343
## Batch2 -0.05599 0.76828 62.00000 -0.073 0.942
## Recov.Days:InitialDomnond -0.02090 0.01932 62.00000 -1.082 0.284
## Recov.Days:Batch2 -0.01323 0.02020 62.00000 -0.655 0.515
## InitialDomnond:Batch2 0.09916 0.94232 62.00000 0.105 0.917
## Recov.Days:InitialDomnond:Batch2 0.01155 0.02403 62.00000 0.480 0.633
##
## Correlation of Fixed Effects:
## (Intr) Rcv.Dy IntlDm Batch2 Rc.D:ID R.D:B2 InD:B2
## Recov.Days -0.796
## InitilDmnnd -0.855 0.680
## Batch2 -0.738 0.587 0.631
## Rcv.Dys:InD 0.706 -0.887 -0.790 -0.521
## Rcv.Dys:Bt2 0.675 -0.848 -0.577 -0.766 0.752
## IntlDmnn:B2 0.602 -0.479 -0.704 -0.815 0.556 0.625
## Rcv.D:ID:B2 -0.568 0.713 0.635 0.644 -0.804 -0.841 -0.757
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular anova(ofavrecovmod2) #no significant interaction between days and batch p=0.515## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Recov.Days 7.5863 7.5863 1 62 5.3699 0.0238 *
## InitialDom 2.1714 2.1714 1 62 1.5370 0.2197
## Batch 0.0003 0.0003 1 62 0.0002 0.9892
## Recov.Days:InitialDom 2.2405 2.2405 1 62 1.5859 0.2126
## Recov.Days:Batch 0.5449 0.5449 1 62 0.3857 0.5368
## InitialDom:Batch 0.0156 0.0156 1 62 0.0111 0.9165
## Recov.Days:InitialDom:Batch 0.3262 0.3262 1 62 0.2309 0.6326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ssidrecovmod<-glmer(Shprop~Recov.Days*InitialDom*Batch+(1|Colony), data=filter(allrecov,allrecov$Species=='S.siderea'))
plot_resqq(ssidrecovmod)
ssidrecovmod2<-as_lmerModLmerTest(ssidrecovmod)
summary(ssidrecovmod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shprop ~ Recov.Days * InitialDom * Batch + (1 | Colony)
## Data: filter(allrecov, allrecov$Species == "S.siderea")
##
## REML criterion at convergence: 403.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8040 -0.1669 -0.0472 0.0603 4.7527
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 2.279 1.51
## Residual 28.516 5.34
## Number of obs: 65, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.11618 2.37121 47.91136 0.049
## Recov.Days 0.01047 0.06825 52.06277 0.153
## InitialDomnond -0.18731 3.07911 47.51902 -0.061
## Batch2 0.40831 2.81832 55.62053 0.145
## Recov.Days:InitialDomnond 0.02058 0.08276 51.92254 0.249
## Recov.Days:Batch2 0.01760 0.07557 52.00907 0.233
## InitialDomnond:Batch2 -10.68440 5.32130 54.93125 -2.008
## Recov.Days:InitialDomnond:Batch2 0.41917 0.11360 52.13169 3.690
## Pr(>|t|)
## (Intercept) 0.961126
## Recov.Days 0.878615
## InitialDomnond 0.951749
## Batch2 0.885332
## Recov.Days:InitialDomnond 0.804567
## Recov.Days:Batch2 0.816757
## InitialDomnond:Batch2 0.049587 *
## Recov.Days:InitialDomnond:Batch2 0.000536 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Rcv.Dy IntlDm Batch2 Rc.D:ID R.D:B2 InD:B2
## Recov.Days -0.748
## InitilDmnnd -0.767 0.576
## Batch2 -0.781 0.632 0.608
## Rcv.Dys:InD 0.617 -0.825 -0.739 -0.520
## Rcv.Dys:Bt2 0.675 -0.903 -0.520 -0.755 0.745
## IntlDmnn:B2 0.413 -0.335 -0.522 -0.532 0.427 0.400
## Rcv.D:ID:B2 -0.449 0.601 0.538 0.502 -0.729 -0.665 -0.763 anova(ssidrecovmod2) #no significant interaction between days and batch p=0.816757## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Recov.Days 639.65 639.65 1 51.848 22.4310 1.729e-05 ***
## InitialDom 108.96 108.96 1 41.211 3.8210 0.0574198 .
## Batch 98.62 98.62 1 54.457 3.4583 0.0683417 .
## Recov.Days:InitialDom 469.32 469.32 1 51.840 16.4579 0.0001680 ***
## Recov.Days:Batch 456.25 456.25 1 52.131 15.9997 0.0002008 ***
## InitialDom:Batch 114.96 114.96 1 54.931 4.0315 0.0495869 *
## Recov.Days:InitialDom:Batch 388.29 388.29 1 52.132 13.6165 0.0005365 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 par(mfrow=c(1,2))
Spec.labs=c('M. cavernosa','O. faveolata','S. siderea')
names(Spec.labs)=c('M.cavernosa','O.faveolata','S.siderea')
all_bleach= ggplot()+
geom_smooth(method='glm',data=filter(allbleach,allbleach$InitialDom=='d',allbleach$Species=='O.faveolata'),aes(y=Shprop,x=DHW,color=..coloursofavD..),method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allbleach,allbleach$InitialDom=='d',allbleach$Species=='S.siderea'),aes(y=Shprop,x=DHW,color=..coloursssidD..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allbleach,allbleach$InitialDom=='nond',allbleach$Species=='S.siderea'),aes(y=Shprop,x=DHW,color=..coloursssidC..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allbleach,allbleach$Species=='M.cavernosa'),aes(y=Shprop,x=DHW,colour=..coloursmcavC..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allbleach,allbleach$InitialDom=='nond',allbleach$Species=='O.faveolata'),aes(y=Shprop,x=DHW,color=..coloursofavB..), method.args = list(family = "quasipoisson"))+
geom_point(data=filter(allbleach, allbleach$PropD!='NA', allbleach$Shprop<1.01), aes(x=DHW, y=Shprop, shape=InitialDom, colour=PropD),position=position_jitter(width=0.5, height=0.1))+
scale_shape_manual(values=c(1,2), guide='none')+
scale_colour_gradient2(low='blue',high='red',mid='purple',midpoint=0.5,breaks=c(0,0.5,1))+
coord_cartesian(clip='off', ylim=c(0,1))+
scale_y_continuous(breaks=c(0,0.5,1))+
theme_minimal(base_size = 13)%+replace%
theme(panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank())+
labs(y='Proportion of Initial Symbiont Density',x='Degree heating weeks',color='Proportion *Durusdinium*:')+
facet_wrap(nrow=3,vars(Species),labeller=labeller(Species=Spec.labs))+
theme(strip.text.x = element_text(face = 'italic'), legend.title=element_markdown(),legend.position = 'top',panel.spacing = unit(1.5, "lines"))
all_bleach
#ggsave('all_bleach.pdf',device='pdf', width=6,height=5)
all_recov= ggplot()+
geom_smooth(method='glm',data=filter(allrecov,allrecov$InitialDom=='d',allrecov$Species=='O.faveolata'),aes(y=Shprop,x=Recov.Days,color=..coloursofavDrecov..),method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allrecov,allrecov$InitialDom=='d',allrecov$Species=='S.siderea'),aes(y=Shprop,x=Recov.Days,color=..coloursssidDrecov..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allrecov,allrecov$InitialDom=='nond',allrecov$Species=='S.siderea'),aes(y=Shprop,x=Recov.Days,color=..coloursssidCrecov..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allrecov,allrecov$Species=='M.cavernosa'),aes(y=Shprop,x=Recov.Days,colour=..coloursmcavCrecov..), method.args = list(family = "quasipoisson"))+
geom_smooth(method='glm',data=filter(allrecov,allrecov$InitialDom=='nond',allrecov$Species=='O.faveolata'),aes(y=Shprop,x=Recov.Days,color=..coloursofavBrecov..), method.args = list(family = "quasipoisson"))+
geom_point(data=filter(allrecov, allrecov$PropD!='NA'), aes(x=Recov.Days, y=Shprop, shape=InitialDom, colour=PropD),position=position_jitter(width=5, height=0.25))+
scale_colour_gradient2(low='blue',high='red',mid='purple', midpoint=0.5,limits=c(0,1),breaks=c(0,0.5,1), guide='none')+
scale_shape_manual(values=c(1,2),labels=c('*Durusdinium*', 'Other'))+
scale_y_continuous(position='right')+
coord_cartesian(clip='on',ylim=c(0,7.5))+
theme_minimal(base_size = 13)%+replace%
theme(panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank())+
labs(y='Proportion of Initial Symbiont Density <br />',x='Days after heat stress', shape='Initially dominant symbiont:')+
facet_wrap(nrow=3,vars(Species),labeller=labeller(Species=Spec.labs))+
theme(axis.title.y.right =element_markdown(),legend.title = element_markdown(),legend.text=element_markdown(),strip.text=element_text(face='italic'),panel.spacing = unit(1.5, "lines"), legend.position = 'top')
all_recov
#ggsave('all_recov.pdf',device='pdf', width=6,height=5)
# UPDATES FOR ROSS: plot now includes recovery day 0 on the right hand side, and a quasipoisson model is fitted. These data are also now relative sh changes. datapoints have also been added. Error bars and 50% gradient colour have both been made more visible. R^2 values are now also reported for these models. QUasipoisson GLMs were plotted in the figure, but for the statistical tests of significance, mixed effects models with poisson distributions did not seem to work well, so I just used gaussian mixed models- is this ok?SHUFFLING
Figure 5: Inter- and intra-species variation in the shuffling towards durusdinium after recovery from heat stress.
Comparison of Proportion D before start of heat stress compared to two months into recovery. Data are binned at intervals of 0.02 (for variable Proportion D), and size aesthetic relates to number of cores in each bin to reduce overplotting. Data from April and October are both included here.
mcavshift<-filter(mcav_wide_batches,!is.na(postDcat),!is.na(preDcat))
ofavshift<-filter(ofav_wide_batches,!is.na(postDcat),!is.na(preDcat))
ssidshift<-filter(ssid_wide_batches,!is.na(postDcat),!is.na(preDcat))
beforeaftershift<-rbind(mcavshift,ofavshift,ssidshift)
beforeaftershift$Treatment<-factor(beforeaftershift$Treatment,levels=c('Manipulated','Control'))
gaindd=expression(paste("Gained",italic(" Durusdinium")))
lostd=expression(paste("Lost",italic(" Durusdinium")))
ggplot()+geom_count(data=(beforeaftershift),
aes(x=preDcat,y=postDcat,colour=Treatment))+
facet_grid(cols=vars(Species),labeller=labeller(Species=Spec.labs))+
geom_abline(slope=1,intercept = 0,linetype='dotted')+
scale_color_manual(values=c('brown2','blue3'),labels=c('Bleached','Control'))+
theme_minimal(base_size = 15)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())+
labs(x='Proportion *Durusdinium* before',y='Proportion *Durusdinium* after')+
guides(size=guide_legend(title='Binned n'))+
scale_size(range=c(1,14),breaks=(c(1,2,6,10)))+
theme(axis.title.x = element_markdown(),axis.title.y=element_markdown())+
theme(strip.text.x = element_text(face = 'italic'),
panel.spacing = unit(1.5, "lines"), legend.position = 'bottom', legend.title=element_text(size=12))+
guides(colour = guide_legend(title=''))
#ggsave('redblue_shift.pdf',device='pdf', width=10,height=5) #add 'gained durusdinium'/'lost durusdinium' labelsFigure 6: Shuffling towards high proportions of Durusdinium after heat stress recovery was greatest in M. cavernosa, followed by O. faveolata, followed by S. siderea.
Emulating the ‘symbiont shuffling’ plot to compare coral species, as in Cunning et al 2018 figure 3b. In order to be able to include mcav, which has no variation in initial proportion d, the following mcav models are independent of initial proportion d, then integrated in with the previous ofav & ssid predicted effects.
shufflemod=lmer(post_propD ~ pre_propD + Treatment*Species*Batch+(1|Colony),
data=batches)
shufflemod2=glm(post_propD~pre_propD + Treatment*Species*Batch, data=batches, family='quasibinomial')
plot_resqq(shufflemod) 
plot(shufflemod2)
#these two models (quasibinomial family vs linear mixed effects) give fairly differnet results...
summary(shufflemod)# Model summary## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_propD ~ pre_propD + Treatment * Species * Batch + (1 | Colony)
## Data: batches
##
## REML criterion at convergence: 38.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70264 -0.41847 -0.03774 0.41929 2.53018
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.01996 0.1413
## Residual 0.05266 0.2295
## Number of obs: 105, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 0.88340 0.07730 53.50301
## pre_propD 0.59718 0.09715 31.47720
## TreatmentControl -0.89499 0.12754 78.92189
## SpeciesO.faveolata -0.25219 0.11305 50.69438
## SpeciesS.siderea -0.42798 0.11937 47.68955
## BatchOctober 0.12362 0.08833 68.59613
## TreatmentControl:SpeciesO.faveolata 0.58869 0.18219 81.18425
## TreatmentControl:SpeciesS.siderea 0.51816 0.18053 77.77864
## TreatmentControl:BatchOctober -0.11775 0.17995 78.02351
## SpeciesO.faveolata:BatchOctober -0.18054 0.13902 72.46015
## SpeciesS.siderea:BatchOctober -0.13107 0.13288 74.16874
## TreatmentControl:SpeciesO.faveolata:BatchOctober -0.10542 0.26833 83.10055
## TreatmentControl:SpeciesS.siderea:BatchOctober 0.36282 0.26333 82.61570
## t value Pr(>|t|)
## (Intercept) 11.428 5.62e-16 ***
## pre_propD 6.147 7.59e-07 ***
## TreatmentControl -7.017 6.95e-10 ***
## SpeciesO.faveolata -2.231 0.03015 *
## SpeciesS.siderea -3.585 0.00079 ***
## BatchOctober 1.400 0.16616
## TreatmentControl:SpeciesO.faveolata 3.231 0.00178 **
## TreatmentControl:SpeciesS.siderea 2.870 0.00528 **
## TreatmentControl:BatchOctober -0.654 0.51482
## SpeciesO.faveolata:BatchOctober -1.299 0.19817
## SpeciesS.siderea:BatchOctober -0.986 0.32718
## TreatmentControl:SpeciesO.faveolata:BatchOctober -0.393 0.69541
## TreatmentControl:SpeciesS.siderea:BatchOctober 1.378 0.17198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#bleaching and initial proportion dusursdinium both had significant effects on shuffling (p=6.95e-10 and p=7.59e-07 respectively). Batch did not significantly affect shuffling (p=0.16616).
#Get fitted values averaging across initial proportion durusdinium
eff <- effect(c('pre_propD', 'Species','Treatment'), shufflemod,
xlevels=list(pre_propD=seq(0, 1, by=0.01)))
# Get all fitted values and subsets for each treatment
res <- droplevels(data.frame(eff))
res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea", "M.cavernosa"))
res.Bl <- subset(data.frame(eff), Treatment=="Manipulated")
res.Ct <- subset(data.frame(eff), Treatment=="Control")
# Get AUC for fitted values, lower and upper confidence limits
auc <- aggregate(res[, c("fit", "lower", "upper")],
by=list(Species=res$Species, Treatment=res$Treatment),
FUN=function(x) (mean(x)-0.5)/0.5)
#force limits of 1 and -1 on ofav and ssid
auc[4,5]=1.0
auc.list <- split(auc, list(auc$Treatment))
auc$Treatment=factor(auc$Treatment, levels=c('Manipulated','Control'))
levels(auc$Treatment)<-c('Bleached','Control')
auc$Species=factor(auc$Species, levels=c('M.cavernosa', 'O.faveolata', 'S.siderea'))
#UPDATES FOR ROSS: I have now changed this so it is based on a mixed effects model (colony as a random factor as in all other analyses), and the result looks fairly differnet from before (which was based on a glm with quasibinomial distribution). Limited predicted values to between 1 and -1. What is the significance of the (-0.5/0.5) in the aggregation function?#using model independent of initial prop d for mcav
shufflemod2=lmer(post_propD ~Treatment*Species*Batch+(1|Colony),
data=batches)
plot_resqq(shufflemod2) 
summary(shufflemod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_propD ~ Treatment * Species * Batch + (1 | Colony)
## Data: batches
##
## REML criterion at convergence: 68.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.56595 -0.29433 0.06883 0.35489 2.28865
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.05966 0.2443
## Residual 0.06123 0.2475
## Number of obs: 108, groups: Colony, 30
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 0.88490 0.10345 44.86844
## TreatmentControl -0.90026 0.14093 75.47542
## SpeciesO.faveolata -0.08228 0.14647 44.85383
## SpeciesS.siderea -0.20658 0.14575 44.55474
## BatchOctober 0.12790 0.09591 69.67963
## TreatmentControl:SpeciesO.faveolata 0.73218 0.20115 77.38500
## TreatmentControl:SpeciesS.siderea 0.57034 0.19799 75.06017
## TreatmentControl:BatchOctober -0.12784 0.19847 74.76893
## SpeciesO.faveolata:BatchOctober -0.14478 0.15171 72.75896
## SpeciesS.siderea:BatchOctober 0.08178 0.13793 71.67140
## TreatmentControl:SpeciesO.faveolata:BatchOctober -0.31940 0.29652 78.95345
## TreatmentControl:SpeciesS.siderea:BatchOctober 0.22211 0.28892 77.88603
## t value Pr(>|t|)
## (Intercept) 8.554 5.59e-11 ***
## TreatmentControl -6.388 1.25e-08 ***
## SpeciesO.faveolata -0.562 0.57707
## SpeciesS.siderea -1.417 0.16334
## BatchOctober 1.334 0.18670
## TreatmentControl:SpeciesO.faveolata 3.640 0.00049 ***
## TreatmentControl:SpeciesS.siderea 2.881 0.00517 **
## TreatmentControl:BatchOctober -0.644 0.52144
## SpeciesO.faveolata:BatchOctober -0.954 0.34308
## SpeciesS.siderea:BatchOctober 0.593 0.55513
## TreatmentControl:SpeciesO.faveolata:BatchOctober -1.077 0.28469
## TreatmentControl:SpeciesS.siderea:BatchOctober 0.769 0.44437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnC SpcsO. SpcsS. BtchOc TrC:SO. TrC:SS. TrC:BO SO.:BO
## TrtmntCntrl -0.371
## SpecsO.fvlt -0.706 0.262
## SpecisS.sdr -0.710 0.263 0.501
## BatchOctobr -0.473 0.396 0.334 0.336
## TrtmntC:SO. 0.260 -0.701 -0.374 -0.185 -0.278
## TrtmntC:SS. 0.264 -0.712 -0.187 -0.362 -0.282 0.499
## TrtmntCn:BO 0.261 -0.718 -0.184 -0.185 -0.538 0.503 0.511
## SpcsO.fv:BO 0.299 -0.250 -0.436 -0.212 -0.632 0.381 0.178 0.340
## SpcsS.sd:BO 0.329 -0.275 -0.232 -0.473 -0.695 0.193 0.408 0.374 0.440
## TrtC:SO.:BO -0.175 0.481 0.265 0.124 0.360 -0.719 -0.342 -0.669 -0.580
## TrtC:SS.:BO -0.179 0.493 0.127 0.264 0.370 -0.346 -0.729 -0.687 -0.234
## SS.:BO TC:SO.:
## TrtmntCntrl
## SpecsO.fvlt
## SpecisS.sdr
## BatchOctobr
## TrtmntC:SO.
## TrtmntC:SS.
## TrtmntCn:BO
## SpcsO.fv:BO
## SpcsS.sd:BO
## TrtC:SO.:BO -0.251
## TrtC:SS.:BO -0.560 0.460#again for the initial durusdinium independent model, batch is not a significant driver of shuffling, p= 0.18670.
eff2 <- effect(c('Species','Treatment','Batch'), shufflemod2)
# Get all fitted values and subsets for each treatment level
res2 <- droplevels(data.frame(eff2))
#force limits of 1 and 0 on mcav confidence intervals
res2[1,7]=1
res2[7,4]=1
res2[7,7]=1
res2[4,4]=0
res2[4,6]=0
res2[10,4]=0
res2[10,6]=0
res2$Species <- factor(res2$Species, levels=c('M.cavernosa',"O.faveolata", "S.siderea"))
res2.Bl <- subset(data.frame(eff2), Treatment=="Manipulated")
res2.Ct <- subset(data.frame(eff2), Treatment=="Control")
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res2[, c("fit", "lower", "upper")],
by=list(Species=res2$Species, Treatment=res2$Treatment),
FUN=function(x) (mean(x)))
levels(auc2$Treatment)<-c('Control','Bleached')
auc2.list <- split(auc2, list(auc2$Treatment))
ofavssidshuff<-auc[c(1,2,4,5),]
mcavshuff<-auc2[c(1,4),]
speciesshuff<-rbind(ofavssidshuff,mcavshuff)
speciesshuff$Species=factor(speciesshuff$Species, levels=c('M.cavernosa', 'O.faveolata', 'S.siderea'))
#quote the shuffling metric for bleached cores of each species, 'S': mcav=0.94245186, ofav=0.80645955, ssid=0.50105511.
ggplot(data=speciesshuff)+
geom_hline(yintercept=0,linetype='dashed')+
geom_errorbar(aes(ymin=lower, ymax=upper, x=Species, colour=Species, group=Treatment),size=0.5, position=position_dodge(width=0.5), width=0.2)+
geom_point(aes(y=fit, x=Species,shape=Treatment, colour=Species),size=2, position=position_dodge(width=0.5))+
theme_minimal(base_size = 13)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank())+
scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
guides(colour=F, shape=guide_legend(title=''))+
theme(legend.position=c(0.1,0.15))+
labs(y='Symbiont Shuffling', x='')+
scale_y_continuous(limits=c(-1,1),expand=c(0,0))+
scale_x_discrete(labels=c(mcav,ofav,ssid), expand=expansion(mult=c(0.5,0.2)))
#ggsave('speciesshuffle.pdf',device='pdf', width=7,height=5)
#add 'more durusdinium/less durusdinium' labels post-save
#QUESTION FOR ROSS: unsure whether or not to use the '-0.5/0.5' in the aggregate function here (it seems to give a more conservative estimate for shuffling in bleached mcav, but does not work for the mcav controls)......as shown below
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res2[, c("fit", "lower", "upper")],
by=list(Species=res2$Species, Treatment=res2$Treatment),
FUN=function(x) (mean(x)-0.5)/0.5)
levels(auc2$Treatment)<-c('Control','Bleached')
auc2## Species Treatment fit lower upper
## 1 M.cavernosa Control -1.0000000000 -1.0000000 -0.4732779
## 2 O.faveolata Control -0.1950366838 -0.7576407 0.3675674
## 3 S.siderea Control 0.0007565134 -0.5582910 0.5598040
## 4 M.cavernosa Bleached 0.8849037164 0.4888545 1.0000000
## 5 O.faveolata Bleached 0.5883544891 0.1428467 1.0338623
## 6 S.siderea Bleached 0.5663219334 0.1568355 0.9758084shuffling and photochemical efficiency
Figure S1: Differnet photochemical disadvantages of hosting Durusdinium at ambient temperatures cannot explain the difference in shuffling between O. faveolata and S. siderea.
Recreating Cunning et al 2018 plot of the relative photochemical disadvantages of hosting durusdinium at ambient temperatures to explain species hierarchy (ofav>ssid) in shuffling.
ipam=filter(allbleach,allbleach$Timepoint=='Pre-bleach',allbleach$Species!='M.cavernosa')
ggplot()+
geom_smooth(method='glm',data=ipam,aes(y=Y2,x=PropD,colour=Species, linetype=Batch),method.args = list(family = "quasibinomial"))+
geom_point(data=ipam,aes(y=Y2,x=PropD,colour=Species, shape=Batch))+
scale_colour_manual(values=c('darkorange1','darkturquoise'),labels=c(ofav,ssid))+
theme_minimal(base_size = 15)+
labs(y='Initial (ambient) Fv/Fm',x='Proportion *Durusdinium*')+
theme(axis.title.x= element_markdown(),legend.title=element_blank(),legend.position='bottom')## `geom_smooth()` using formula 'y ~ x'
#ggsave('inity2.pdf',device='pdf', width=7,height=5)
ipaminitmod=lmer(Y2~PropD*Species*Batch+(1|Colony),data=ipam)
plot_resqq(ipaminitmod)
summary(ipaminitmod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y2 ~ PropD * Species * Batch + (1 | Colony)
## Data: ipam
##
## REML criterion at convergence: -203.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.86465 -0.49506 -0.03011 0.50441 1.80649
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0001145 0.01070
## Residual 0.0002989 0.01729
## Number of obs: 51, groups: Colony, 19
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.5293324 0.0070538 25.0437420 75.042 < 2e-16
## PropD -0.0032823 0.0134075 27.9645585 -0.245 0.808391
## SpeciesS.siderea 0.0001896 0.0108844 26.2466965 0.017 0.986235
## Batch2 0.0363517 0.0096538 37.8622808 3.766 0.000564
## PropD:SpeciesS.siderea -0.0216401 0.0186729 31.1020055 -1.159 0.255317
## PropD:Batch2 -0.0017877 0.0163661 33.5505728 -0.109 0.913669
## SpeciesS.siderea:Batch2 -0.0099099 0.0158472 33.9985376 -0.625 0.535923
## PropD:SpeciesS.siderea:Batch2 -0.0053723 0.0239426 32.1793643 -0.224 0.823879
##
## (Intercept) ***
## PropD
## SpeciesS.siderea
## Batch2 ***
## PropD:SpeciesS.siderea
## PropD:Batch2
## SpeciesS.siderea:Batch2
## PropD:SpeciesS.siderea:Batch2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PropD SpcsS. Batch2 PrD:SS. PrD:B2 SS.:B2
## PropD -0.554
## SpecisS.sdr -0.648 0.359
## Batch2 -0.517 0.289 0.335
## PrpD:SpcsS. 0.398 -0.718 -0.621 -0.207
## PropD:Btch2 0.317 -0.564 -0.205 -0.634 0.405
## SpcsS.sd:B2 0.315 -0.176 -0.420 -0.609 0.246 0.386
## PrpD:SS.:B2 -0.217 0.386 0.312 0.434 -0.526 -0.684 -0.731#UPDATE FOR ROSS: I'm not sure this plot is needed, maybe just the stats, given there isn't a significant finding here. No significant differnece in the interaction between proportion durusdinium and coral species on Fv/Fm (p=0.255). However, initial Fv/Fm values were generally significantly higher in October (compared to april), regardless of coral species and proportion durusdinium (p=0.0005), although I don't think this is particularly relevant to the story of this paper?Figure S2: Despite greater symbiont loss and reduced symbiont recovery in the October M. cavernosa corals (figure 4), shuffling between April and October were similar for all three coral species.
#before getting fitted values, test statistical significance of batch for each species inidivually:
ofavbatchmod=lmer(post_propD ~ pre_propD + Batch+(1|Colony),
data=filter(batches, batches$Species=='O.faveolata', Treatment=='Manipulated'))
ssidbatchmod=lmer(post_propD ~ pre_propD + Batch+(1|Colony),
data=filter(batches, batches$Species=='S.siderea', Treatment=='Manipulated'))
mcavbatchmod=lmer(post_propD ~Batch+(1|Colony),
data=filter(batches, batches$Species=='M.cavernosa', Treatment=='Manipulated'))
summary(ofavbatchmod) #p=0.546701## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_propD ~ pre_propD + Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "O.faveolata", Treatment ==
## "Manipulated")
##
## REML criterion at convergence: 11.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.84265 -0.10735 0.09364 0.20476 1.67041
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.08313 0.2883
## Residual 0.03990 0.1997
## Number of obs: 22, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.70531 0.12833 7.45285 5.496 0.000736 ***
## pre_propD 0.35253 0.23111 6.33618 1.525 0.175414
## BatchOctober -0.06274 0.10098 11.32362 -0.621 0.546701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pr_prD
## pre_propD -0.539
## BatchOctobr -0.202 -0.145summary(ssidbatchmod) #p=0.85918## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_propD ~ pre_propD + Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "S.siderea", Treatment ==
## "Manipulated")
##
## REML criterion at convergence: 17.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15335 -0.08276 -0.01649 0.13550 1.81479
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.03293 0.1815
## Residual 0.07492 0.2737
## Number of obs: 25, groups: Colony, 9
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.50616 0.13036 11.20440 3.883 0.00247 **
## pre_propD 0.48309 0.19318 11.98865 2.501 0.02790 *
## BatchOctober 0.02277 0.12658 18.68964 0.180 0.85918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pr_prD
## pre_propD -0.658
## BatchOctobr -0.101 -0.420summary(mcavbatchmod) #p=0.0505## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: post_propD ~ Batch + (1 | Colony)
## Data: filter(batches, batches$Species == "M.cavernosa", Treatment ==
## "Manipulated")
##
## REML criterion at convergence: -10.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7890 -0.2055 0.1068 0.4819 1.2017
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.01003 0.1002
## Residual 0.02478 0.1574
## Number of obs: 28, groups: Colony, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.88288 0.05391 18.20147 16.377 2.43e-12 ***
## BatchOctober 0.12636 0.06082 20.57959 2.078 0.0505 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BatchOctobr -0.576eff <- effect(c('pre_propD', 'Species','Treatment', 'Batch'), shufflemod,
xlevels=list(pre_propD=seq(0, 1, by=0.01)))
# Get all fitted values and subsets for each treatment and now also batch
res <- droplevels(data.frame(eff))
res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea", "M.cavernosa"))
res.Bl <- subset(data.frame(eff), Treatment=="Manipulated")
res.Ct <- subset(data.frame(eff), Treatment=="Control")
res.Ap<- subset(data.frame(eff), Batch=='April')
res.Oc<- subset(data.frame(eff), Batch=='October')
# Get AUC for fitted values, lower and upper confidence limits
auc <- aggregate(res[, c("fit", "lower", "upper")],
by=list(Species=res$Species, Treatment=res$Treatment, Batch=res$Batch),
FUN=function(x) (mean(x)-0.5)/0.5)
auc[4,6]=1.0
auc[10,6]=1.0
auc$Treatment=factor(auc$Treatment, levels=c('Manipulated','Control'))
levels(auc$Treatment)<-c('Bleached','Control')
auc$Species=factor(auc$Species, levels=c('M.cavernosa', 'O.faveolata', 'S.siderea'))
auc$Batch=factor(auc$Batch, levels=c('April','October'))
ofavssidbatchshuff<- auc[-c(3,6,9,12),]
eff2 <- effect(c('Species','Treatment','Batch'), shufflemod2)## NOTE: SpeciesTreatmentBatch is not a high-order term in the modelres2 <- droplevels(data.frame(eff2))
#force limits of 1 and 0 on mcav confidence intervals
res2[1,7]=1
res2[7,4]=1
res2[7,7]=1
res2[4,4]=0
res2[4,6]=0
res2[10,4]=0
res2[10,6]=0
res2$Species <- factor(res2$Species, levels=c('M.cavernosa',"O.faveolata", "S.siderea"))
res2.Bl <- subset(data.frame(eff2), Treatment=="Manipulated")
res2.Ct <- subset(data.frame(eff2), Treatment=="Control")
res2.Ap<- subset(data.frame(eff2), Batch=='April')
res2.Oc<- subset(data.frame(eff2), Batch=='October')
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res2[, c("fit", "lower", "upper")],
by=list(Species=res2$Species, Treatment=res2$Treatment, Batch=res2$Batch),
FUN=function(x) (mean(x)))
levels(auc2$Treatment)<-c('Control','Bleached')
mcavbatchshuff<- auc2[c(1,4,7,10),]
batchshuff<- rbind(ofavssidbatchshuff, mcavbatchshuff)
batchshuff$spbatch<- interaction (batchshuff$Batch,batchshuff$Species)
ggplot(data=filter(batchshuff,batchshuff$Treatment=='Bleached'))+
geom_errorbar(aes(ymin=lower, ymax=upper, x=Species, colour=Species, group=Batch),size=0.5, position=position_dodge(width=0.5), width=0.2)+
geom_point(aes(y=fit, x=Species, shape=Batch, colour=Species),size=2, position=position_dodge(width=0.5))+
theme_minimal(base_size = 13)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank())+
scale_color_manual(values=c('deeppink2','darkorange1', 'darkturquoise'))+
guides(colour=F, shape=guide_legend(title=''))+
theme(legend.position=c(0.1,0.15))+
labs(y='Symbiont Shuffling in Bleached Corals', x='')+
scale_y_continuous(limits=c(0,1),expand=c(0,0))
#ggsave('shufflebatches.pdf',device='pdf',width=7,height=5)Figure 7: Shuffling occurred later in recovery in M. cavernosa compared to O. faveolata and S. siderea.
shuffletimes=read.csv('shuffle_timepoints.csv',header =T)
shuffletimes$Treatment=factor(shuffletimes$Treatment)
shuff=read.csv('shuff.csv',header = T)
shuff$Colony=factor(shuff$Colony)
shuff$Species=factor(shuff$Species)
shuff$Timepoint=as.integer(shuff$Timepoint,length=3)
shuffmod=lmer(PropD~PropD0+Species*Timepoint+(1|Colony),data=shuff)
plot_resqq(shuffmod)
summary(shuffmod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PropD ~ PropD0 + Species * Timepoint + (1 | Colony)
## Data: shuff
##
## REML criterion at convergence: 104.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5126 -0.4678 0.0141 0.4194 3.7174
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.03418 0.1849
## Residual 0.07168 0.2677
## Number of obs: 208, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.42076 0.09910 92.91601 -4.246 5.16e-05
## PropD0 0.39576 0.09866 45.80985 4.011 0.000221
## SpeciesO.faveolata 0.77968 0.14591 96.40165 5.344 6.09e-07
## SpeciesS.siderea 0.94733 0.16291 97.27006 5.815 7.71e-08
## Timepoint 0.44927 0.03626 174.74719 12.391 < 2e-16
## SpeciesO.faveolata:Timepoint -0.32474 0.05375 174.54374 -6.042 8.96e-09
## SpeciesS.siderea:Timepoint -0.42711 0.05572 175.65745 -7.665 1.17e-12
##
## (Intercept) ***
## PropD0 ***
## SpeciesO.faveolata ***
## SpeciesS.siderea ***
## Timepoint ***
## SpeciesO.faveolata:Timepoint ***
## SpeciesS.siderea:Timepoint ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PropD0 SpcsO. SpcsS. Timpnt SpO.:T
## PropD0 -0.002
## SpecsO.fvlt -0.679 -0.202
## SpecisS.sdr -0.608 -0.358 0.485
## Timepoint -0.736 0.000 0.500 0.447
## SpcsO.fvl:T 0.496 -0.013 -0.729 -0.297 -0.675
## SpcsS.sdr:T 0.479 0.012 -0.327 -0.720 -0.651 0.439anova(shuffmod,test='F')## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PropD0 1.1534 1.1534 1 45.81 16.091 0.0002209 ***
## Species 3.0202 1.5101 2 100.17 21.066 2.308e-08 ***
## Timepoint 5.4417 5.4417 1 175.23 75.914 2.180e-15 ***
## Species:Timepoint 4.8481 2.4241 2 175.19 33.817 3.808e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1propD1=lmer(PropD ~ PropD0+Timepoint*Species+(1|Colony),
data=shuff)
eff1 <- effect(c('PropD0', 'Timepoint', 'Species'), propD1,
xlevels=list(PropD0=seq(0, 1, by=0.01)))
res <- droplevels(data.frame(eff1))
res$Species <- factor(res$Species, levels=c("O.faveolata", "S.siderea"))
# Get AUC for fitted values, lower and upper confidence limits
auc1 <- aggregate(res[, c("fit", "lower", "upper")],
by=list(Species=res$Species,Timepoint=res$Timepoint),
FUN=function(x) (mean(x)-0.5)/0.5)
#now an intial prop d independent model for mcav
propD2=lmer(PropD ~ Timepoint*Species+(1|Colony),
data=shuff)
eff2 <- effect(c('Timepoint', 'Species'), propD2)## NOTE: TimepointSpecies is not a high-order term in the modelres <- droplevels(data.frame(eff2))
res$Species <- factor(res$Species, levels=c("M.cavernosa"))
# Get AUC for fitted values, lower and upper confidence limits
auc2 <- aggregate(res[, c("fit", "lower", "upper")],
by=list(Species=res$Species,Timepoint=res$Timepoint),
FUN=function(x) (mean(x)))
aucall <- rbind(auc1, auc2)
aucall$Timepoint=factor(aucall$Timepoint, levels=c(1,2,3))
aucall[11,4]=0
aucall[9,5]=1
aucall[15,5]=1
aucall=aucall[-c(3,4,7,8,12,14),]
aucall$Species=factor(aucall$Species, levels=c('M.cavernosa','O.faveolata','S.siderea'))
##########
ofav=expression(paste(italic("O. faveolata")))
ssid=expression(paste(italic("S. siderea")))
mcav=expression(paste(italic("M. cavernosa")))
ggplot(aucall, aes(x = Timepoint, y = fit, group = Species)) +
geom_errorbar(data=aucall, aes(ymin = lower, ymax = upper),
position = position_dodge(0.2), lwd = 0.2, width = 0.2) +
geom_point(aes(color = Species),
position = position_dodge(0.2), size = 2.5)+
geom_hline(yintercept = 0, lwd = 0.1) +
scale_y_continuous(limits = c(-0.25, 1))+
scale_x_discrete(labels=c('Post heat stress','1 month recovery','2 month recovery'),expand=expansion(mult=c(0.2,0.4)))+
theme_minimal(base_size = 13)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank())+
geom_line(linetype='dashed',alpha=0.6,aes(colour=Species),position = position_dodge(0.2))+
scale_colour_manual(values=c('deeppink2','darkorange1','darkturquoise'),labels=c(mcav,ofav,ssid), name='')+
labs(y='Cumulative Symbiont Shuffling',x='')+
annotate(geom='text',x=3.4,y=0.1,label=gaindd,size=4)+
annotate(geom='text',x=3.4,y=-0.1,label=lostd,size=4)+
theme(legend.position = 'bottom')
#ggsave('shuffletimingcumulative.pdf',device='pdf',height=5,width=7)
###UPDATE FOR ROSS: again, I have changed this to usea mixed effects model, and have forced limits of 1 and -1 for the errorbars (1 and 0 for mcav which had no initial durusdinium). The predicted values are based of a model which controls for pre-heat stress proportion durusdinium (rather than proportion durusdinium at the previous timepoint). Again, I'm unsure of the '-0.5/0.5' part of the aggregate function- this plot shows lower predicted shuffling in ofav and ssid compared to original plots based on a glm with quasibinomial distribution. Making 'timepoint' an integer rather than a factor in the model also changed things. Would liek to perform stats on the interaction effect on shuffling between timepoint and species, but unsire how since mcav fitted to a separate model. Figure 8: O. faveolata and S.siderea corals that initially hosted background Durusdinium shuffled more quickly and more fully to Durusdinium than those with no initial Durusdinium.
#now look at timing for ofav and ssid switching and shuffling
shufftimesofss=read.csv('shuffle_timepointsofavssid.csv',head=T)
shufftimesofss=filter(shufftimesofss,(shufftimesofss$Treatment=='Manipulated'))
shufftimesofss=filter(shufftimesofss,(shufftimesofss$initial.d!='NA'))
shuffmod2=lmer(PropD2~initial.d*Species+(1|Colony),data=shufftimesofss)
plot_resqq(shuffmod2)
summary(shuffmod2) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PropD2 ~ initial.d * Species + (1 | Colony)
## Data: shufftimesofss
##
## REML criterion at convergence: 14.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06511 -0.05183 -0.03832 0.37651 1.61534
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.07282 0.2699
## Residual 0.06080 0.2466
## Number of obs: 19, groups: Colony, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.04340 0.15703 7.73680 0.276 0.78952
## initial.dy 0.59597 0.19799 14.11762 3.010 0.00929 **
## SpeciesS.siderea -0.04340 0.30244 9.04668 -0.143 0.88905
## initial.dy:SpeciesS.siderea -0.09597 0.41571 10.65192 -0.231 0.82180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intl.d SpcsS.
## initial.dy -0.633
## SpecisS.sdr -0.519 0.329
## intl.dy:SS. 0.302 -0.476 -0.688anova(shuffmod2,test='F')# prop d significantly higher in those that started with some d p=0.02371## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## initial.d 0.42255 0.42255 1 10.6519 6.9504 0.02371 *
## Species 0.01054 0.01054 1 8.0401 0.1733 0.68804
## initial.d:Species 0.00324 0.00324 1 10.6519 0.0533 0.82180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shuffmod3=lmer(PropD3~initial.d*Species+(1|Colony),data=shufftimesofss)
plot_resqq(shuffmod3)
summary(shuffmod3) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PropD3 ~ initial.d * Species + (1 | Colony)
## Data: shufftimesofss
##
## REML criterion at convergence: 24.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2491 -0.8478 0.0000 0.9146 1.2343
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0000 0.0000
## Residual 0.1641 0.4051
## Number of obs: 22, groups: Colony, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.506007 0.143226 18.000000 3.533 0.00238 **
## initial.dy 0.493993 0.202552 18.000000 2.439 0.02532 *
## SpeciesS.siderea -0.006007 0.248074 18.000000 -0.024 0.98095
## initial.dy:SpeciesS.siderea 0.006007 0.405104 18.000000 0.015 0.98833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intl.d SpcsS.
## initial.dy -0.707
## SpecisS.sdr -0.577 0.408
## intl.dy:SS. 0.354 -0.500 -0.612
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularanova(shuffmod3,test='F')# prop d significantly higher in those that started with some d p=0.02456## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## initial.d 0.98802 0.98802 1 18 6.0205 0.02456 *
## Species 0.00004 0.00004 1 18 0.0002 0.98833
## initial.d:Species 0.00004 0.00004 1 18 0.0002 0.98833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shuffmod4=lmer(PropD4~initial.d*Species+(1|Colony),data=shufftimesofss)
plot_resqq(shuffmod4)
summary(shuffmod4)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PropD4 ~ initial.d * Species + (1 | Colony)
## Data: shufftimesofss
##
## REML criterion at convergence: 20.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.22228 -0.70611 -0.04998 0.15747 1.72820
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02019 0.1421
## Residual 0.11494 0.3390
## Number of obs: 22, groups: Colony, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.48417 0.15660 10.12194 3.092 0.0113 *
## initial.dy 0.46078 0.20573 11.85230 2.240 0.0451 *
## SpeciesS.siderea -0.10025 0.22646 7.52952 -0.443 0.6704
## initial.dy:SpeciesS.siderea 0.08997 0.36966 12.24422 0.243 0.8117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intl.d SpcsS.
## initial.dy -0.731
## SpecisS.sdr -0.692 0.506
## intl.dy:SS. 0.407 -0.557 -0.601anova(shuffmod4,test='F') # prop d significantly higher in those that started with some d p=0.01777## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## initial.d 0.86064 0.86064 1 12.244 7.4878 0.01777 *
## Species 0.00999 0.00999 1 9.893 0.0869 0.77421
## initial.d:Species 0.00681 0.00681 1 12.244 0.0592 0.81173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1shufftimesofss= pivot_longer(data=shufftimesofss, cols=starts_with('PropD'),
names_to='Timepoint', names_prefix='PropD', values_to='PropD',
values_drop_na=T)
switchex=expression(paste('no initial ', italic("Durusdinium ")))
shuffex=expression(paste('>0 initial ', italic("Durusdinium "), '≤ 25%'))
ggplot(shufftimesofss, aes(x = Timepoint, y = PropD, group=initial.d)) +
stat_summary(aes(shape=initial.d, colour=initial.d),fun.data='mean_se',
position = position_dodge(0.2), size=0.5) +
scale_x_discrete(labels=c('Pre heat stress','Post heat stress','1 month recovery','2 month recovery'))+
scale_shape_manual(values=c(0,15),labels=c(switchex,shuffex))+
scale_colour_manual(values=c('blue3','brown2'),labels=c(switchex,shuffex))+
stat_summary(geom='line', aes(linetype=initial.d, colour=initial.d),position = position_dodge(0.2))+
scale_linetype_manual(values=c('dashed','solid'),labels=c(switchex,shuffex))+
theme_minimal(base_size = 15)%+replace%
theme(panel.border = element_rect(colour='grey20',fill=NA),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank())+
labs(y="Proportion *Durusdinium* <br /> (*O. faveolata* & *S. siderea*)", x='')+
theme(legend.position = 'bottom',legend.title=element_text(size=0), axis.title.y = element_markdown())
#ggsave('ofavshuffswitch.pdf',device='pdf',height=5,width=7)
###UPDATE FOR ROSS: I have corrected this analysis to be based on proportion d (mean +-SE), rather than the shuffling metric, given that this is a comparison between those with and without any initial d. Consistent with figure 4, 25% initial d was taken as the threshold, so those with no initial d were compared against those with >0 but <25% initial d. Due to small sample size, ofav and ssid are grouped together here (model finds no significant affect of species). Other questions for Ross: 1. Light intensity calculation. 2. Based on a reviewer comment replacement of D. trenchii with Durusdinium.